/*
 *  simulation.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *
 */

#include "simulation.h"

Simulation::Simulation()
:	nProc(1),
	nSys(0),
	ini(),
//	fPos(),
	keyword(),
	inputFileName(),
	title(),
	dimension(3),
	mode(0),
	vecType(0),
	matType(0),
	vProcesses(),
	mSysProc(),
	vAnalyses(),
	currentAnalysis(),
	mMaterials(),
	mSurfaces(),
	rawtime(),
	vDofNames(),
	vStressNames(),
	graphicsFormat(),
	graphicsPath("./"),
	graphicsStep(0.0),
	fGraphTot(),
	fEner(),
	fEntropy(),
	printPosFlag(0),
	printEnerFlag(0),
	printIterFlag(0),
	printTempFlag(0),
	printEntropyFlag(0),
	printAMFlag(0),
	printConFlag(0),
	printPProjFlag(0),
	printStressFlag(0),
	plotDispFlag(0),
	plotProcFlag(0),
	plotStressFlag(0),
	outEnerFile(),
	outEntropyFile(),
	outTempFile(),
	Ekin(0.0),
	Epot(0.0),
	Econ(0.0),
	Eamb(0.0),
	Dint(0.0),
	Entropy(0.0),
	EntropyAmb(0.0),
	AM(),
	vCon(),
	vTempNames(),
	symmFlag(true)
{}

// ***

Simulation::Simulation(char* fileName)
:	nProc(1),
	nSys(0),
	ini(),
//	fPos(),
	keyword(),
	inputFileName(),
	title(),
	dimension(3),
	mode(0),
	vecType(0),
	matType(0),
	vProcesses(),
	mSysProc(),
	vAnalyses(),
	currentAnalysis(),
	mMaterials(),
	mSurfaces(),
	rawtime(),
	vDofNames(),
	vStressNames(),
	graphicsFormat(),
	graphicsPath("./"),
	graphicsStep(0.0),
	fGraphTot(),
	fEner(),
	printPosFlag(0),
	printEnerFlag(0),
	printIterFlag(0),
	printTempFlag(0),
	printEntropyFlag(0),
	printAMFlag(0),
	printConFlag(0),
	printPProjFlag(0),
	printStressFlag(0),
	plotDispFlag(0),
	plotProcFlag(0),
	plotStressFlag(0),
	outEnerFile(),
	outTempFile(),
	Ekin(0.0),
	Epot(0.0),
	Econ(0.0),
	Eamb(0.0),
	Dint(0.0),
	AM(),
	vCon(),
	vTempNames(),
	symmFlag(true)
{
	AM.resize(dimension);
	vCon.resize(3);
	
    strcpy(inputFileName,fileName);
	
	readGeneral();		/** Purpose: read general options from input */ 
	
	if(!keyword.compare("BEGINSYSTEM")) readSystems();
						/** Purpose: Read systems from input file */
						/** */
						/** 1) Creates processes */
						/** 2) Within each process, creates systems */
						/** 3) Within each system, creates a list of its points (with coordinates), but */
						/**	   it does not create the actual points (objects) */
						/** 4) Creates bodies. At each body, creates components. */
						/** At the end of this stage, each component has a list with the names of its points */
						/** Some of these names may belong to the local process, others may be remote */
	
	if(!keyword.compare("BEGINMATERIAL")) readMaterials();	/** Read materials and send them to processes */

	if(!keyword.compare("BEGINSURFACE")) readSurfaces();	/** Read surfaces and send them to processes */

	if(!keyword.compare("BEGINPOST")) readPost();			/** Read postprocess options */
	
	if(!keyword.compare("BEGINANALYSIS")) readAnalyses();
						/** General purpose: Read analyses from input file and sequentially perform calculations */
						/** */
						/** 1) Creates an object analysis */
						/** 2) Create/remove loads */
						/** 3) Create/remove systems/bodies/components/points */
						/** 4) Send the analysis to each process. Each process reconstruct the analysis object  */
						/** 5) For any analysis that adds/remove some component/body/system (and for the very first analysis) */
						/**	   perform the following steps: */
						/**		- Reset system tree. TO BE DONE
						/**		- Build a tree of systems within each process (creating a root system) */
						/**		- Set pointers to materials at components */
						/**		- Create the points (objects) within each system */
						/**		- At each component, set pointers to their actual points, based on its list of point names. */
						/**		  Going through this list of point names:
						/**			(i) If the name belongs to a local point, sets the pointer to it */
						/**			(ii) If the name belongs to a remote point:	+ create a point at the root system; */
						/**														+ set the pointer to this (now local) point */
						/**														+ add its name in a list located at the remote process */
						/**		- Set the crossed relations among processes. */
						/**		   At this stage, each process knows to which remote processes has to communicate with,*/
						/**		   and what are the shared points. */
						/**		   This is done by means of: */
						/**			(i) Creation of interface objects within each process. */
						/**			    Each interface relates two processes, and has a list of names of shared points (local and nonlocal) */
						/**			(ii) Each interface gets a pointer to its related remote interface */
						/**			(iii) Related interfaces are synchronized: shared point names are the same and in the same order */
						/**		- At each body, set the names and pointers of its points */
						/**		- Set dofs of points, with the info given by components */
						/**		- At each point, set global numbers for its dofs */
						/**		- At each component, set the local-global numbering of dofs */
						/**	6) If the analysis is not the first, and does not add/remove any component: */
						/**		- If loads have been added: */
						/**			i) set pointers to points at loads */
						/**			ii)set local/global dof numbering at components / loads */
						/**		- If not, this means that only a change of analysis happened. Check if velocity is involved in the new one */
						/**	7) Perform calculation (loop on t if dynamics) */
						/** 8) Get solution (global dofs)	*/
	double tend = MPI::Wtime();
	lProcess->writeToFileTime("Total execution time : ", tend-(lProcess->tini));
}

// ***

Simulation::~Simulation() { ini.close(); }

// ***

void Simulation::readGeneral()
{

	time(&rawtime);
	ini.open(inputFileName);
	
	while(ini >> keyword && keyword.compare("BEGINSYSTEM") && keyword.compare("BEGINMATERIAL") && keyword.compare("BEGINANALYSIS")  ) 
	{
		if(!keyword.compare("TITLE"))
		{
			ini >> title ;
		}
		else if(!keyword.compare("DIMENSION"))
		{
			ini >> dimension ;
			if (dimension <=0 || dimension > 3)
			{
				cout << "*ERROR: Invalid space dimension \n";
			}
		}
		else if(!keyword.compare("MODE"))
		{
			ini >> mode ;
			// mode=0	File out-processNumber.dat is created with the main results of the calculation
			// mode=1	No files are created
			// mode=2	Just files with time information are created
		}
		else if(!keyword.compare("VECTORTYPE"))
		{
			ini >> vecType ;
		}
		else if(!keyword.compare("MATRIXTYPE"))
		{
			ini >> matType ;
		}
		else if(!keyword.compare("PROCESSES"))
		{
			ini >> nProc ;
			/** Process creation */
			if(nProc>=1)
			{
				for(int k=0; k<nProc; k++)	/** Process numbers are 0,1,... */
				{
					vProcesses.push_back( (new Process(k, dimension, vecType, matType, mode, title)) );
				}
			} else cout << "*ERROR: At least one process must be defined \n";
		}
	}
	if(rankk==0)
	{
		cout << "---------------------------------- \n";	
		cout << "READING GENERAL SIMULATION OPTIONS \n";	
		cout << "---------------------------------- \n";
	}
	cout << "Mode = " << mode << endl;
	if(mode==0)
	{
		// Write general options to output file
		int i;
		for(i=0; i<nProc; i++)
		{
			vProcesses[i]->writeToFile("EPPI OUTPUT FOR PROCESS ",i);
			vProcesses[i]->writeToFile("----------------------------");
			vProcesses[i]->writeToFile("The current local time is: ", ctime (&rawtime));
			vProcesses[i]->writeToFile("----------------------------");
			vProcesses[i]->writeToFile("GENERAL OPTIONS ");
			vProcesses[i]->writeToFile("----------------");
			vProcesses[i]->writeToFile("TITLE: ",title);
			vProcesses[i]->writeToFile("DIMENSION: ",dimension);
			vProcesses[i]->writeToFile("MODE: ",mode);
			vProcesses[i]->writeToFile("NUMBER OF PROCESSES: ",nProc);
			vProcesses[i]->writeToFile("-------------------------");
		}
	}
}

// ***

void Simulation::readSystems()
{
	string auxs;					/** Auxiliary string */
	int auxi;						/** Auxiliary integer */
    pair<string, int> mSysProcEntry;

	if(rankk==0)
	{
		cout << "----------------\n";	
		cout << "READING SYSTEMS \n";	
		cout << "----------------\n";
	}
	
	if(mode==0)
	{
		for(auxi=0; auxi<nProc; auxi++)
		{
			vProcesses[auxi]->writeToFile("READING SYSTEMS");
			vProcesses[auxi]->writeToFile("----------------");
		}
	}

	while(ini >> keyword &&  keyword.compare("BEGINMATERIAL") && keyword.compare("BEGINANALYSIS") && keyword.compare("BEGINSURFACE")  ) 
	{
		if(!keyword.compare("PROC"))
		{
			ini >> auxi;
			if(auxi>=0 & auxi<nProc)
			{
				ini >> auxs ;				/** Read the system name */
				if(!auxs.compare("NAME"))
				{
					ini >> auxs;
				}	else cout << "*ERROR: System defined without name \n";
				if( vProcesses[auxi]->readSystem(auxs,ini,dimension) )	/** Success in telling the process to read the system from the input stream */
				{
					mSysProcEntry.first = auxs;
					mSysProcEntry.second = auxi;
					mSysProc.insert(mSysProcEntry);						
				} else cout << "*ERROR: System " << auxs <<  " has not been added to process " << auxi << "\n" << "----------\n";
			}
			else cout << "*ERROR: Process number " << auxi << " does not exist, there are " << vProcesses.size() << " processes defined\n";
		}
		sendSysProc();		/** Current global system-process table is sent to each process */
	}
	nSys = mSysProc.size();
}

// ***

void Simulation::readAnalyses()
{
	string aName="";					/** Auxiliary string */
	string tName="";					/** Auxiliary string */
	double tEnd=0.0;
	int i;

	if(rankk==0)
	{
		cout << "---------------- \n";	
		cout << "READING ANALYSES \n";	
		cout << "---------------- \n";
	}

	if(mode==0)
	{
		for(i=0; i<nProc; i++)
		{
			vProcesses[i]->writeToFile("--------------------------------");
			vProcesses[i]->writeToFile("READING ANALYSIS");
			vProcesses[i]->writeToFile("--------------------------------");
		}
	}

	while(ini >> keyword && ini.eof()==false )
	{
		if(!keyword.compare("NAME"))
		{
			ini >> aName;
		}
		else if(!keyword.compare("TYPE"))
		{
			ini >> tName;
		}
		if(aName.compare("") && tName.compare(""))
		{
			if(rankk==0) cout << "Reading data from Analysis " << aName << ", type " << tName << " ...\n";
			if(mode==0)
			{
				for(i=0; i<vProcesses.size(); i++)	vProcesses[i]->writeToFile(aName,", type ",tName);
			}

			if(!tName.compare("static")) vAnalyses.push_back(new Static(ini,aName,tName,dimension,graphicsFormat,graphicsPath, graphicsStep,this)) ;
			if(!tName.compare("directdynamic")) vAnalyses.push_back(new DirectDynamic(ini,aName,tName,dimension,tEnd,graphicsFormat,graphicsPath, graphicsStep,  this)) ;
			currentAnalysis = vAnalyses.back();
			if(currentAnalysis==vAnalyses.front())	/** If this is the first analysis ... */
			{
				currentAnalysis->setRebuild(true);		/** Rebuild flag true */
				currentAnalysis->setChangeLoads(true);	/** ChangeLoads flag true  */
			}
			initializeAnalysis();		/** Initialize current analysis */
			tEnd = currentAnalysis->calculate();			/** Perform calculation. Returns the end time (0.0 if not time-dependent) */
			
			if (rankk==0) cout << "----------\n";
			if (rankk==0) cout << "Analysis " << aName << " performed.\n";
			if (rankk==0) cout << "----------\n";
			aName="";
			tName="";
//			if(vDofNames.empty()==false)	writeFileDofs();
			// if(graphicsFormat.compare(""))	writeFileGraph(); //edit: graphic file is now written in static()
		}
	}
	
	if(graphicsFormat.compare("") && rankk==0)	createTotalGraphOutFile();	// Just for the first process
}

// ***

void Simulation::initializeAnalysis()
{
	int i;
	
	if( currentAnalysis->haveToRebuild==true )
	{
		initializeProc(); /** If the processes have to be completely rebuild */
		eraseProcesses();	/** erase processes that do not match rankk.  Only lProcess in left */
		lProcess-> setSizeInterfaces(currentAnalysis->implicitFlag);
		setOutput();
		// if (graphicsFormat.compare("")) createGraphOutFile(); // edit: graphic file is now in static()
	}
	else if ( currentAnalysis->changeLoads==true )
	{
		for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersLoadBodies();				/** Set pointers to points at loads in bodies */
//		for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersLoadJoints();				/** Set pointers to points at loads in joints */
		for(i=0; i<nProc; i++)	vProcesses[i]->setLocalGlob();		/** Set local-global numbering on components */
		checkVelNeeded();	
	}
	
//	else checkVelNeeded();			/** Check if velocity vector dq is needed (for consecutive analyses) */
}

// ***

void Simulation::checkVelNeeded()
{	
//	int i;
//	for(i=0; i<vProcesses.size(); i++)	vProcesses[i].checkVelNeeded();
}

// ***

void Simulation::sendMaterials()
{
	for(int j=0; j<nProc; j++)	vProcesses[j]->setMaterials(mMaterials);
}

// ***

void Simulation::sendSysProc()
{
	int i;
	vector<Process*> vp(nProc);
	for(i=0; i<nProc; i++)	vp[i] = (vProcesses[i]);

	for(i=0; i<nProc; i++)	vProcesses[i]->setSysProc(mSysProc,vp);	/** Send tables to processes */		
}

// ***

void Simulation::readMaterials()
{
	int i;
	string matName="";
	string materialType="";
	pair<string, Material*> matEntry;
	string auxs;
	Material* auxMat;

	if(rankk==0)
	{
		cout << "---------------------------------------------- \n";	
		cout << "READING MATERIALS								\n";	
		cout << "---------------------------------------------- \n";
	}
	
	if(mode==0)
	{
		for(i=0; i<nProc; i++)
		{
			vProcesses[i]->writeToFile("--------------------------------");
			vProcesses[i]->writeToFile("READING MATERIALS");
			vProcesses[i]->writeToFile("--------------------------------");
		}
	}
	
	while(ini >> keyword && keyword.compare("BEGINPOST")!=0 && keyword.compare("BEGINSURFACE") ) 
	{
		if(!keyword.compare("NAME"))
		{
			ini >> matName;
		}
		else if(!keyword.compare("TYPE"))
		{
			ini >> materialType;
		}
		if(matName.compare("") && materialType.compare(""))
		{
			if (rankk==0) cout << "Reading data from Material " << matName << ", type " << materialType << " ...\n";
			if(mode==0)
			{
				for(i=0; i<nProc; i++)
				{
					vProcesses[i]->writeToFile(matName, ", type ", materialType);
				}
			}
			matEntry.first = matName;
			if(!materialType.compare("elasticlog")) matEntry.second = new ElasticLog(ini,matName,materialType);
			else if(!materialType.compare("elasticlinear")) matEntry.second = new ElasticLinear(ini,matName,materialType);
			else if(!materialType.compare("svk")) matEntry.second = new Svk(ini,matName,materialType);
			else if(!materialType.compare("neohooke")) matEntry.second = new NeoHooke(ini,matName,materialType);
            else if(!materialType.compare("mooneyrivlin")) matEntry.second = new MooneyRivlin(ini,matName,materialType);
			else if(!materialType.compare("elasticlog3d")) matEntry.second = new ElasticLog3D(ini,matName,materialType);
			else if(!materialType.compare("ogden")) matEntry.second = new Ogden(ini,matName,materialType);
            else if(!materialType.compare("composite"))
			{
				// Read the name of the elastic material, in order to send the pointer
				ini >> auxs;	// Reads the ELASTICMAT label
				ini >> auxs;	// name of the elastic material
				auxMat = mMaterials[auxs];
				matEntry.second = new  Composite(ini,matName,materialType,auxMat);
			}
            else if(!materialType.compare("damage"))
			{
				// Read the name of the elastic material, in order to send the pointer
				ini >> auxs;	// Reads the ELASTICMAT label
				ini >> auxs;	// name of the elastic material
				auxMat = mMaterials[auxs];
				matEntry.second = new  Damage(ini,matName,materialType,auxMat);
			}
			else if(!materialType.compare("viscoelastic")) 
			{
				// Read the name of the elastic material, in order to send the pointer
				ini >> auxs;	// Reads the ELASTICMAT label
				ini >> auxs;	// name of the elastic material
				auxMat = mMaterials[auxs];
				matEntry.second = new Viscoelastic(ini,matName,materialType,auxMat);
			}
			else if(!materialType.compare("thermoelastic"))
			{
				// Read the name of the elastic material, in order to send the pointer
				ini >> auxs;	// Reads the ELASTICMAT label
				ini >> auxs;	// name of the elastic material
				auxMat = mMaterials[auxs];
				matEntry.second = new  Thermoelastic(ini,matName,materialType,auxMat);
			}
			mMaterials.insert(matEntry);
			matName="";
			materialType="";

		}
	}
	sendMaterials();
}

// ***

void Simulation::readPost()
{
    int i;
	int nPostPos;
	int nTemp;
	string pointName;
	string bodyName;
	string line;
        
	stringstream stnum;

	if(rankk==0)
	{
		cout << "--------------------------- \n";	
		cout << "READING POSTPROCESS OPTIONS \n";	
		cout << "--------------------------- \n";
	}

	if(mode==0)
	{
		for(i=0; i<nProc; i++)
		{	
			vProcesses[i]->writeToFile("--------------------------------");
			vProcesses[i]->writeToFile("READING POSTPROCESS OPTIONS");
			vProcesses[i]->writeToFile("--------------------------------");
		}
	}
	while(ini >> keyword && keyword.compare("BEGINANALYSIS")!=0)
	{
		if(!keyword.compare("PRINTDOFS"))
		{
//			ini >> nPostPos;
			
			getline(ini,line);
			//cout << "line is " << line << "\n";
			istringstream iss(line);
			while (iss >> keyword )
			{
				pointName = keyword;
				vDofNames.push_back(pointName);
				//cout << "keyword is " << keyword << endl;
			}
			
			
			/*
			for (i=0; i<nPostPos; ++i)
			{
				ini >> pointName;
				vDofNames.push_back(pointName);
			}
			*/
			printPosFlag = 1;
		}
		if(!keyword.compare("PRINTSTRESS"))
		{
			getline(ini,line);
			//cout << "line is " << line << "\n";
			istringstream iss(line);
			while (iss >> keyword )
			{
				pointName = keyword;
				//cout << "keyword is " << keyword << endl;
				vStressNames.push_back(pointName);
			}
			printStressFlag = 1;			
		}
		if(!keyword.compare("GRAPHICS"))
		{
			std::vector<string> gFormatParam;

			getline(ini,line);			// Reads the type of 3D tool
			istringstream iss(line);

			while (iss >> keyword )
			{
				gFormatParam.push_back(keyword);
			}

			graphicsFormat = gFormatParam[0];
			if(gFormatParam.size()>1)
			{
				stnum << gFormatParam[1];
				stnum >> graphicsStep;
			}

			//cout << "graphicsFormat=" << graphicsFormat << endl;
			//cout << "graphicsStep=" << graphicsStep << endl;
		}
		if(!keyword.compare("GRAPHICSPATH"))
		{
			ini >> graphicsPath;
                        
			if( !checkForSlash(graphicsPath))
			{
			      char *currentdir = new char[300];

			      currentdir = getcwd(currentdir,300); // get current directory's path

			      graphicsPath = std::string(currentdir) + "/" + graphicsPath + "/"; // within current directory's path
			      //graphicsPath = "./" + graphicsPath + "/"; // within current directory's path
			}

			struct stat st;
			if(stat(graphicsPath.c_str(),&st) != 0)
				mkdir(graphicsPath.c_str(), S_IRWXU);
                        
			cout << "Graphics data will be saved in " << graphicsPath << endl;
		}
		if(!keyword.compare("PRINTENER"))
		{
			printEnerFlag = 1;
			createFileEner();
		}
		if(!keyword.compare("PRINTITER"))
		{
			printIterFlag = 1;
		}		
		if(!keyword.compare("PRINTTEMPERATURE"))
		{
			getline(ini,line);
			istringstream iss(line);
			while (iss >> keyword )
			{
				bodyName = keyword;
				vTempNames.push_back(bodyName);
			}
			printTempFlag = 1;						
		}
		if(!keyword.compare("PRINTENTROPY"))
		{
			printEntropyFlag = 1;
			createFileEntropy();
		}		
		if(!keyword.compare("PRINTAMOMENTUM"))
		{
			printAMFlag = 1;
			createFileAM();
		}
		if(!keyword.compare("PRINTCONSTRAINT"))
		{
			printConFlag = 1;
			createFileCon();
		}
		if(!keyword.compare("PRINTPROJPENALTY"))
		{
			printPProjFlag = 1;
			createFilePProj();
		}
		if(!keyword.compare("PLOTDISP"))
		{
			plotDispFlag = 1;
		}
		if(!keyword.compare("PLOTPROC"))
		{
			plotProcFlag = 1;
		}
		if(!keyword.compare("PLOTSTRESS"))
		{
			plotStressFlag = 1;
		}				
	}
}

// ***

void Simulation::initializeProc()
{
	int i;
	cout << "Initializing processes ...\n";

	// FIRST, IT IS NECESSARY TO RESET SOME OBJECTS WITHIN EACH PROCESS, SYSTEM,...  TO BE DONE */
	// 
	//**

	for(i=0; i<nProc; i++)	vProcesses[i]->buildSysTree();	/** Build system tree */
	for(i=0; i<nProc; i++)	vProcesses[i]->setMatAtComp();	/** Set pointer to material at components */
	for(i=0; i<nProc; i++)	vProcesses[i]->setSurfAtBodies();	/** Set pointer to surface at bodies */
//	for(i=0; i<nProc; i++)	vProcesses[i]->setLocalPointsAtBodies();	/** Set local points at bodies (just coordinates, these points do not have dofs) */
	for(i=0; i<nProc; i++)	vProcesses[i]->repairChangedPointNames();
	for(i=0; i<nProc; i++)	vProcesses[i]->createContacts();
//	for(i=0; i<nProc; i++)	vProcesses[i]->setBodiesAtJoints();
	for(i=0; i<nProc; i++)	vProcesses[i]->createPoints();	/** Create points within systems */

	if (rankk==0) cout << "	Setting crossed relations among processes ...\n";
	for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersCompBodies();	/** Set pointers to points at components in bodies */
	for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersLoadBodies();	/** Set pointers to points at loads in bodies */
//	for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersCompJoints();	/** Set pointers to points at components in joints */
//	for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersLoadJoints();	/** Set pointers to points at loads in joints */
	for(i=0; i<nProc; i++)	vProcesses[i]->setPPointersConsJoints();	/** Set pointers to points at constraints in joints */
	for(i=0; i<nProc; i++)	vProcesses[i]->setPPointSurface();			/** Set pointers to surface points in bodies */
	for(i=0; i<nProc; i++)	vProcesses[i]->setRelProc();			/** Set process relations */
	for(i=0; i<nProc; i++)	vProcesses[i]->createInterfaces(currentAnalysis->integratorType);	/** Creates interfaces */
	for(i=0; i<nProc; i++)	vProcesses[i]->relateInterfaces();	/** Setup and synchronize interfaces */
	
	if (rankk==0) cout << "	Setting body points lists and looking for orphan points ...\n";	/** NEEDS REVISON, only detects locally to each process */
	if (rankk==0) cout << "	Setting degrees of freedom of processes ...\n";

	/** Set the dynFlag at processes */
	if( currentAnalysis->type.compare("static") )	for(i=0; i<nProc; i++)	vProcesses[i]->dynFlag=1;
	else for(i=0; i<nProc; i++)	vProcesses[i]->dynFlag=0;
	// Set the projFlag at processes
	if (currentAnalysis->projFlag==true)	for(i=0; i<nProc; i++)	vProcesses[i]->projFlag=true;
	
	for(i=0; i<nProc; i++)	vProcesses[i]->setDofsBodies();		/** Set dofs inside points in bodies, from info given by components */
	for(i=0; i<nProc; i++)	vProcesses[i]->setDofsRoots();		/** Set dofs inside points at root systems. */
	for(i=0; i<nProc; i++)	vProcesses[i]->setDofsJoints();		/** Set dofs inside points in joints, from info given by components */
	
//	cout << "A-0" << endl;
	if (rankk==0) cout << "	Setting global dof numbers ...\n";
	for(i=0; i<nProc; i++)	vProcesses[i]->setGlobalDofNumbers();	/** Count and set global dof numbers */
//	cout << "A-1" << endl;
	if (rankk==0) cout << "	Setting local-global numbering ...\n";
	for(i=0; i<nProc; i++)	vProcesses[i]->setLocalGlob();			/** Set local-global numbering on components */
//	cout << "A-2" << endl;
	for(i=0; i<nProc; i++)	vProcesses[i]->calcSizeInterfaces(currentAnalysis->integratorType, currentAnalysis->implicitFlag);	/** Give proper dof size to interfaces and set local/global numberings */
	if(mode==0) reportContentFile();
	/** Set ambient temperature */
	for(i=0; i<nProc; i++)	vProcesses[i]->setTemperature();
	/**  Initialization of the dofs of bodies at each process */
	for(i=0; i<nProc; i++)	vProcesses[i]->initializeGlobalDofsBodies();
	// Initialize dofs that have not been initialized at bodies
	// This could be the case of joints relating local and remote points.
	for(i=0; i<nProc; i++)	vProcesses[i]->synchronizeSharedDofs();
	/**  Initialization of the dofs of joints at each process */
	for(i=0; i<nProc; i++)	vProcesses[i]->initializeGlobalDofsJoints();
	// Set the vector containing the imposed dof's displacements
	for(i=0; i<nProc; i++)	vProcesses[i]->setNoDofs();
	for(i=0; i<nProc; i++)	vProcesses[i]->setImposedDofs();
	for(i=0; i<nProc; i++)	vProcesses[i]->setHistories(mHistories);	
	
	/** Set initial conditions */
	for(i=0; i<nProc; i++)	vProcesses[i]->setInitialConditions(currentAnalysis->h);
	/** Set q (and dq if required) at components */
	for(i=0; i<nProc; i++)	vProcesses[i]->setQatComponents();
	
	
	if (currentAnalysis->lumpedMassFlag==true)
	{
		// Create lumped mass matrices (for explicit integration)
		for(i=0; i<nProc; i++)	vProcesses[i]->createLumpedMass();	
		// Synchronize lumped mass matrices (for explicit integration)
		for(i=0; i<nProc; i++)	vProcesses[i]->synchronizeLumpedMass();	
	}
	
	cout << "Initialization of processes completed ...\n";
}
	
// ***

void Simulation::reportContentFile()
{
	int i;
	for(i=0; i<nProc; i++)
	{
		vProcesses[i]->writeToFile("-------");
		vProcesses[i]->writeToFile("SUMMARY");
		vProcesses[i]->writeToFile("-------");
		vProcesses[i]->writeToFile("Number of DOFs : ",vProcesses[i]->numDofs);
		vProcesses[i]->writeToFile("Point Name \t DOFs");
		vProcesses[i]->writeToFile("-----------   -------------------");
		vProcesses[i]->reportPointsDofsFile();
		vProcesses[i]->writeToFile("--------------");
		vProcesses[i]->writeToFile("END OF SUMMARY");
		vProcesses[i]->writeToFile("--------------");
	}
}

// ***

void Simulation::eraseProcesses()
{
	if(vProcesses.size()>=1)
	{
		lProcess = vProcesses[rankk];
		vProcesses.erase(vProcesses.begin(),vProcesses.end());
		vector<Process*>(vProcesses).swap(vProcesses);
	}	
}

// ***

void Simulation::setOutput()
{
	int i;
	int numProc=0;

	// DOF output
	if (printPosFlag==1)
	{
		for(i=0; i<vDofNames.size(); ++i)	// Loop on points
		{
			numProc = mSysProc[getStringHead(vDofNames[i])];
			if( numProc==rankk)	// Identify the process where the point belongs
			{
				lProcess->createDofOutFile();
				lProcess->fillVDofsOut( vDofNames[i] );
			}
		}
		currentAnalysis->printDofsFlag = 1;
	}
	// STRESS output
	if (printStressFlag==1)
	{
		for(i=0; i<vStressNames.size(); ++i)	// Loop on points
		{
			numProc = mSysProc[getStringHead(vStressNames[i])];
			if( numProc==rankk)	// Identify the process where the point belongs
			{
				lProcess->createStressOutFile();
				lProcess->fillVStressOut( vStressNames[i] );
			}
		}
		currentAnalysis->printStressFlag = 1;
	}	
	// ENERGY output
	if (printEnerFlag ==1)
	{
		currentAnalysis->printEnerFlag = 1;
	}
	// ITERATIONS output
	if (printIterFlag==1)
	{
		currentAnalysis->printIterFlag = 1;
	}
	// TEMPERATURE output
	if (printTempFlag ==1)
	{
		for(i=0; i<vTempNames.size(); ++i)	// Loop on bodies/joints
		{
			numProc = mSysProc[getStringHead(vTempNames[i])];
			if( numProc==rankk)	// Identify the process where the body belongs
			{
				lProcess->createTempOutFile();
				lProcess->fillVTempOut( vTempNames[i] );
			}
		}		
		currentAnalysis->printTempFlag = 1;
	}	
	// ENTROPY output
	if (printEntropyFlag ==1)
	{
		currentAnalysis->printEntropyFlag = 1;
	}
	// ANGULAR MOMENTUM output
	if (printAMFlag ==1)
	{
		currentAnalysis->printAMFlag = 1;
	}
	// CONSTRAINT output
	if (printConFlag ==1)
	{
		currentAnalysis->printConFlag = 1;
	}
	// PROJECTION PENALTY output
	if (printPProjFlag ==1)
	{
		currentAnalysis->printPProjFlag = 1;
	}
	// DISPLACEMENT plot
	if (plotDispFlag ==1)
	{
		currentAnalysis->plotDispFlag = 1;
	}
	// PROCESS plot
	if (plotProcFlag ==1)
	{
		currentAnalysis->plotProcFlag = 1;
	}
	// STRESS plot
	if (plotStressFlag ==1)
	{
		currentAnalysis->plotStressFlag = 1;
	}
	// TIME output
	if (mode == 2 || mode == 0)
	{
		currentAnalysis->printTimeFlag = 1;
	}						
	
}

// ***

void Simulation::setGravity(double& grav, int& gravDir)
{
	int i;
	for(i=0; i<nProc; i++)	vProcesses[i]->setGravity(grav, gravDir);
}

// ***
void Simulation::createTotalGraphOutFile()
{
	char outGraphFile[300];
	char outGraphFileP[300];
	
	if(!graphicsFormat.compare("vtu"))	
	{
		cout << "Creating total vtu output file ...\n";
		int i;
		string stp;
		stringstream stnumproc;
		string snumproc;
		ifstream fgp;
		string auxs;

		strcpy(outGraphFile, graphicsPath.c_str());
		strcat(outGraphFile,"GraphTot.pvd");

		fGraphTot.open(outGraphFile);
		fGraphTot << "<?xml version=\"1.0\"?>\n";
		fGraphTot << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
		fGraphTot << "	<Collection>\n";

		// Open pvd files
		strcpy(outGraphFileP,graphicsPath.c_str());
		for(i=0; i<nProc;++i)
		{
			cout << "Writting results of process " << i << "...\n";
			stnumproc << i;
			stnumproc >> snumproc;

			stp = "Graph"+snumproc+".pvd";
			strcat(outGraphFileP,stp.c_str());
			fgp.open(outGraphFileP);
			getline(fgp,auxs);
			getline(fgp,auxs);
			getline(fgp,auxs);
			while(auxs.compare("	</Collection>"))
			{
				getline(fgp,auxs);
				if (auxs.compare("	</Collection>")) fGraphTot << auxs << "\n";
			}
			stnumproc.clear();
			snumproc.clear();
			fgp.close();
		}
		fGraphTot << "	</Collection>\n";
		fGraphTot << "</VTKFile>";
		fGraphTot.close();
		cout << "Graphical output file finished \n";
	}	
}

// ***

void Simulation::collectEner()
{
	Ekin = 0.0;
	Epot = 0.0;
	Econ = 0.0;
    Eext = 0.0;
	Eamb = 0.0;
    Dint = 0.0;
        
	MPI::COMM_WORLD.Allreduce(&(lProcess->Ekin),&Ekin,1,MPI::DOUBLE,MPI::SUM);
	MPI::COMM_WORLD.Allreduce(&(lProcess->Epot),&Epot,1,MPI::DOUBLE,MPI::SUM);
	MPI::COMM_WORLD.Allreduce(&(lProcess->Econ),&Econ,1,MPI::DOUBLE,MPI::SUM);
    MPI::COMM_WORLD.Allreduce(&(lProcess->Eext),&Eext,1,MPI::DOUBLE,MPI::SUM);
	MPI::COMM_WORLD.Allreduce(&(lProcess->Eamb),&Eamb,1,MPI::DOUBLE,MPI::SUM);
    MPI::COMM_WORLD.Allreduce(&(lProcess->Dint),&Dint,1,MPI::DOUBLE,MPI::SUM);
}

// ***

void Simulation::collectEntropy()
{
	Entropy = 0.0;
	EntropyAmb = 0.0;
	
	MPI::COMM_WORLD.Allreduce(&(lProcess->Entropy),&Entropy,1,MPI::DOUBLE,MPI::SUM);
	MPI::COMM_WORLD.Allreduce(&(lProcess->EntropyAmb),&EntropyAmb,1,MPI::DOUBLE,MPI::SUM);
}

// ***

void Simulation::collectAM()
{
	setZero(AM);
	
	MPI::COMM_WORLD.Allreduce(&(lProcess->AM[0]),&AM[0],dimension,MPI::DOUBLE,MPI::SUM);
}

// ***

void Simulation::collectCon()
{
	setZero(vCon);
	MPI::COMM_WORLD.Allreduce(&(lProcess->vCon[0]),&vCon[0],3,MPI::DOUBLE,MPI::SUM);
}

// ***

void Simulation::createFileEner()
{	
	// Create if the output file does not exist
	if ( !strcmp(outEnerFile,"") )
	{
		strcpy(outEnerFile, graphicsPath.c_str());
		strcat(outEnerFile,"Ener.dat");

		//cout << outEnerFile << endl;

		fEner.open(outEnerFile);
	}	
}

// ***

void Simulation::createFileEntropy()
{	
		strcpy(outEntropyFile,graphicsPath.c_str());
		strcat(outEntropyFile,"Entropy.dat");
		
		fEntropy.open(outEntropyFile);
}

// ***

void Simulation::createFileAM()
{	
	strcpy(outAMFile,graphicsPath.c_str());
	strcat(outAMFile,"AMomentum.dat");
	fAM.open(outAMFile);
}

// ***

void Simulation::createFileCon()
{	
	strcat(outConFile,graphicsPath.c_str());
	strcpy(outConFile,"Constraint.dat");
	fCon.open(outConFile);
}

// ***

void Simulation::createFilePProj()
{	
	strcat(outPProjFile,graphicsPath.c_str());
	strcpy(outPProjFile,"PProj.dat");
	fPProj.open(outPProjFile);
}

// ***

void Simulation::writeFileEner(double& t)
{
	fEner << t <<" "<< Ekin <<" "<< Epot <<" "<< Econ <<" "<< Ekin+Epot <<" "<< Ekin+Epot+Econ <<" "<< Eamb <<" "<< Ekin+Epot+Econ+Eamb <<" "<< Eext <<" "<< Ekin+Epot+Econ+Eext <<" "<< Dint << endl;
}

// ***

void Simulation::writeFileEntropy(double& t)
{
	fEntropy << t << " " << Entropy << " " << EntropyAmb <<  " " << Entropy+EntropyAmb << endl;
}

// ***

void Simulation::writeFileAM(double& t)
{
	int i;
	int mod = 0.0;
	
	fAM << t;
	for(i=0; i<dimension; ++i)
	{
		fAM << " " << AM[i];
		mod+= AM[i]*AM[i];
	}
	fAM << " " << sqrt(mod);
	fAM << endl;
}

// ***

void Simulation::writeFileCon(double& t)
{
	int i;
	int mod = 0.0;
	
	fCon << t;
	
	// Square root of squares
	for(i=0; i<3; ++i)
	{
		fCon << " " << sqrt(vCon[i]);
	}
	fCon << endl;
}

// ***

void Simulation::writeFilePProj(double& t, double& pproj)
{
	fPProj << t;
	fPProj << " " << pproj << endl;
}

// ***

void Simulation::closeFileEner(void)
{
	fEner.close();
}

// ***

void Simulation::closeFileEntropy(void)
{
	fEntropy.close();
}

// ***

void Simulation::closeFileAM(void)
{
	fAM.close();
}

// ***

void Simulation::closeFileCon(void)
{
	fCon.close();
}

// ***

void Simulation::closeFilePProj(void)
{
	fPProj.close();
}

// ***

void Simulation::setProcTemperature(double& theta)
{
	int i;
	for(i=0; i<nProc; i++)	vProcesses[i]->thetab = theta;
}

// ***

void Simulation::readSurfaces()
{
	int i;
	string surfName="";
	string surfType="";
	
	if(rankk==0)
	{
		cout << "---------------------------------------------- \n";	
		cout << "READING SURFACES								\n";	
		cout << "---------------------------------------------- \n";
	}
	
	if(mode==0)
	{
		for(i=0; i<nProc; i++)
		{
			vProcesses[i]->writeToFile("--------------------------------");
			vProcesses[i]->writeToFile("READING SURFACES");
			vProcesses[i]->writeToFile("--------------------------------");
		}
	}
	
	while(ini >> keyword && keyword.compare("BEGINPOST") ) 
	{
		if(!keyword.compare("NAME"))
		{
			ini >> surfName;
		}
		else if(!keyword.compare("TYPE"))
		{
			ini >> surfType;
		}
		if(surfName.compare("") && surfType.compare(""))
		{
			if (rankk==0) cout << "Reading data from Surface " << surfName << ", type " << surfType << " ...\n";
			if(mode==0)
			{
				for(i=0; i<nProc; i++)
				{
					vProcesses[i]->writeToFile(surfName, ", type ", surfType);
				}
			}
			if (!surfType.compare("mesh"))
			{
				mSurfaces[surfName] = new SurfaceMesh(ini,surfName);
			}
			else if (!surfType.compare("sphere"))
			{
				mSurfaces[surfName] = new SurfaceSphere(ini,surfName);
			}
			else if (!surfType.compare("ellipsoid"))
			{
				mSurfaces[surfName] = new SurfaceEllipsoid(ini,surfName);
			}
			else if (!surfType.compare("superellipsoid"))
			{
				mSurfaces[surfName] = new SurfaceSuperEllipsoid(ini,surfName);
			}
			else if (!surfType.compare("plane"))
			{
				mSurfaces[surfName] = new SurfacePlane(ini,surfName);
			}			
			else if (!surfType.compare("box"))
			{
				mSurfaces[surfName] = new SurfaceBox(ini,surfName);
			}			
			
			surfName="";
			surfType="";
		}
	}
	sendSurfaces();
}

// ***

void Simulation::sendSurfaces()
{
	for(int j=0; j<nProc; j++)	vProcesses[j]->setSurfaces(mSurfaces);
}

// ***

void Simulation::createHistory(char* fName)
{
	string fileName = fName;
	mHistories.insert(pair<string,History*>(fileName,new History(fName)));
}

// ***

void Simulation::createHistory(double& value)
{
	string fileName = "value";
	stringstream stnum;
	string snum;
	stnum << value;
	stnum >> snum;
	fileName+="-";
	fileName+=snum;
	
	//cout << "filename is " << fileName << endl;
	
	mHistories.insert(pair<string,History*>(fileName,new History(value)));
}

